---
title: "Replication File for Minilateralism and Backlash at Security Studies"
author: "Leah Matchett"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

library(dplyr)
library(lubridate)
library(ggplot2)
library(ggthemes)
library(countrycode)
library(Matching)
library(fastDummies)
library(Rcpp)
library(MatchIt)
library(Zelig)
library(tidyr)
library(stargazer)
library(cobalt)
library(rgenoud)
library(extrafont)
font_import()

setwd('Data')
nss = read.csv("NSSintermediate.csv")
#covars = read.csv('NSS_countries_covars.csv')
covars = read.csv('NSS_countries_covars_full.csv')
unideal = read.csv('IdealpointestimatesAll_Apr2020.csv')
cinc = read.csv('NMC_5_0.csv')
gdp = read.csv('API_NY.GDP.MKTP.KD_DS2_en_csv_v2_1865143.csv')

```

#Database Prep

##Import covariates, create needed variables
```{r}
#_______________________________________________________________
#___________ UN ideal points____________________________________
#_______________________________________________________________
# Converting session votes to year. Technically some sessions are early the next year because votes roll over but Im going to go with 1945=1
unideal$year = 1945+ unideal$session

# calculate average UN ideal points for the 2005-2009 period (5 years before the NSS)
un_subset= unideal[unideal$year>2004 & unideal$year<2010,]
unmeans = un_subset %>% group_by(ccode) %>% summarize(avg = mean(IdealPoint, na.rm=TRUE) )
unmeans$idealptavg = unmeans$avg #renaming for the merge


# add countrycodes to NSS intermediate
ccodes = countrycode(nss$country, origin = "country.name", destination = "cown")
nss$ccode = ccodes
# Serbia is not matched - goes by Yugoslavia (345)
nss$ccode[nss$country=="Serbia"]= 345

#merge in UNideal points by year
test = merge(nss, unideal[,c("ccode","year", "IdealPoint")], by.x=c("ccode","year"), by.y = c("ccode","year"), all.x=TRUE)
nss = test

#merge in avg unideal points
test = merge(nss, unmeans[,c("ccode","idealptavg")], by="ccode",all.x = TRUE)
nss = test

#_______________________________________________________________
#___________________ GDP data from WB___________________________
#_______________________________________________________________
# looking only at 2016 year
#gdp$Country.Name = gdp$ï..Country.Name # This can sometimes be needed on windows

gdp$ccode = countrycode(gdp$Country.Name, origin = "country.name",destination = "cown")
# most of the errors are the WB regions, can be ignored. Only adding in the one that matters for my analysis
gdp$ccode[gdp$Country.Name=="Serbia"] = 345

#merge to intermediate df
test = merge(nss, gdp[,c('ccode','GDP')], by="ccode",all.x=TRUE)
nss= test

#______________________________________________________________
#--------------CINC for population data _______________________
#______________________________________________________________

cinc2 = cinc[cinc$year ==2012,] # subset to most recent year to avoid duplicates

test = merge(nss, cinc2[,c('ccode','tpop')], by = c('ccode'), all.x=TRUE)
nss=test

#______________________________________________________________
#_____________Treaty membership data __________________________
#______________________________________________________________

#covars$country = covars$ï..country # This can sometimes be needed on windows
# add ccodes
covars$ccode = countrycode(covars$country, origin = "country.name",destination="cown")
# Errors S√£o Tom√© and Pr√≠ncipe, Serbia : manually fixing
covars$ccode[covars$country=="Serbia"]= 345
covars$ccode[covars$country=="S√£o Tom√© and Pr√≠ncipe"]= 403

#merge to nss df
test= merge(nss, covars, by="ccode", all.x=TRUE)
nss = test

nss$year= nss$year -2011 # need to reset these so 2012 = 1, 2013= 2, etc

nss2 = nss[,c("Region.y","NWS","NPP","IdealPoint","idealptavg","tpop","GDP","NAM","BanTreaty","CTBTrat","PSI","GTRI","GICNT","GP","ICSANT","CPPNM2005","NWFZ","ATT","NATO","NumSupplyGp","NSS2012.2016.x","globalnorms", "country.x", "year", "ccode")]

nss2= nss2[complete.cases(nss2),]

#______________________________________________________________
#_____________Other Clearning __________________________
#______________________________________________________________

#Make region a factor
nss2$Region.y = nss2$Region.y %>% as.factor()
nss3 = nss2

#nss3 = nss2[nss2$NWS<2,] # These are such outliers and all included, it makes more sense to exclude them

# Log Population and GDP
nss3$pop = log(nss3$tpop) 
nss3$GDP= log(nss3$GDP)

nss3$LtnAm = 0 # Make reference Level
nss3 = dummy_cols(nss3, select_columns = "Region.y")



### Add the 2012 NTI global norms value by country, for appendix robustness check
norms2012 = nss3[nss3$year==1,]
norms2012$gn2012 = norms2012$globalnorms 

test = merge(nss3, norms2012[,c("ccode","gn2012")], by="ccode", all.x=T )
nss3= test


#______________________________________________________________
# _____________________Add policy similarity to US_____________
#______________________________________________________________

# add distance from US in UN idealpoint
nss3$USdist= 2.74569- nss3$idealptavg # higher values will be farther away


#cleanup
rm(unideal, unmeans,un_subset, covars, cinc, cinc2, gdp, test)
```



# Section 3 Rhetoric:

##Additional Variables and Rhetorical matching 
```{r}
#--- add the rhetorical DV

iaea = read.csv('Data/IAEAStatements.csv')
#iaea$country = iaea$ï..country #needed for windows sometimes

iaea$ccode = countrycode(iaea$country, origin="country.name",destination = "cown")
iaea$ccode[iaea$country=="Serbia"]= 345

test = merge(nss3, iaea[,c("ccode","SeptStatement","DecStatement")], by = "ccode", all.x=TRUE)
nss3 = test

#--- collapse into one year (2016)
nss_2016 = nss3[nss3$year==5,]

# first need to remove all countries that don't have statements
nss2016_s = nss_2016[!is.na(nss_2016$SeptStatement)|!is.na(nss_2016$DecStatement) ,] # needs to have a value in one column

# Combined Statement Measure: there are a limited number of statements available, so I have combined these into a score of the total positive/negative statement towards the NSS in both September and December Conferences. E.g. a score of 2 would be two positive statements, a score of -2 would be two negative statements. 

nss2016_s$statement = rowSums(nss2016_s[,c("SeptStatement","DecStatement")], na.rm=T)
nss2016_s$statement # 104 observations
nss2016_s$NSS2012.2016.x %>% summary()# 48 % treated

mout3 = matchit(NSS2012.2016.x~ NWS+NPP+idealptavg+pop+ GDP+ NAM + CTBTrat+ PSI + GTRI+ GICNT+ GP+ CPPNM2005+ NWFZ+ NumSupplyGp + Region.y_1 + Region.y_2 + Region.y_3 + Region.y_4 + Region.y_5 + Region.y_6 + Region.y_7 + Region.y_8,  data = nss2016_s, method= "genetic", pop.size = 5000, replace=TRUE)

bal.tab(mout3)#balance table

```

## Rhetoric Stargazer Model
```{r}

rhet_mod = zelig(statement~ NSS2012.2016.x + NWS+NPP+ USdist+ GDP+ NAM + CTBTrat+ PSI + GTRI+ GICNT+ GP+ CPPNM2005+ NWFZ+ NumSupplyGp + NATO+ Region.y_1 + Region.y_2 + Region.y_3 + Region.y_4 + Region.y_5 + Region.y_6 + Region.y_7+ NATO, data = match.data(mout3), model = "ls", robust=T, cluster="ccode")


### Table 3
stargazer( lm(rhet_mod$zelig.out$z.out[[1]]), 
           dep.var.labels=c("Rhetorical Treatment of NSS"),
           covariate.labels=c("Included at NSS","Nuclear Weapons State","Nuclear Power","Dist. from US","GDP","Non-aligned Movement","CTBT ratification","PSI","GTRI","GICNT","G8GP","CPPNM 2005 Amendment","NWFZ","Number of Supply Groups", "NATO","Africa","Asia","Europe-EU","Europe-NonEU","Latin America","North America","MENA","Oceana"), align=TRUE,title= "Zelig Model Regression Results: Rhetoric", nospace= TRUE , single.row = TRUE, star.cutoffs = c(0.05, 0.01, 0.001), notes = c(" * p<0.05; ** p<0.01; *** p<0.001"), notes.append = F, digits = 2)


```


## Figure 2:  Rhetoric Box and Whisker Plot
```{r}

# Simulation ----  Appendix plot A2.4 Figure 3
x.out3 <- setx(rhet_mod, NSS2012.2016.x=0)
x1.out3 <- setx(rhet_mod, NSS2012.2016.x=1)
s.out3 <- sim(rhet_mod, x = x.out3, x1 = x1.out3)

par(mar=c(1,1,1,1))

pdf(file="Output/zeligFigure3.pdf")
plot(s.out3)# Appendix plot A2.4 Figure 3
dev.off()


summary(s.out3) #Data is manually imported from this to the below box plot

mat = matrix(NA, 2, 2)

df = as.data.frame(mat)

min(nss_2016$statement[nss$NSS2012.2016.x==0], na.rm=TRUE)# min unincluded
max(nss_2016$statement[nss$NSS2012.2016.x==0], na.rm=TRUE)

min(nss_2016$statement[nss$NSS2012.2016.x==1], na.rm=TRUE)# min unincluded
max(nss_2016$statement[nss$NSS2012.2016.x==1], na.rm=TRUE)

# these numbers are taken from the min and max variables calculated above
df$ymin = c(-2,-2)
df$ymax = c(2,2)

# these numbers are taken directly from summary(s.out3)
df$low = c(-1.009874, 0.684442)
df$mid = c(-0.2052314, 0.9932506 )
df$top = c(0.5486814, 1.30306)

df$included = c("Excluded","Included")
df$included = factor(df$included, levels = c("Excluded","Included"))

library(extrafont)

font_import()
#font_import( paths = '/Users/leah/Library/Fonts') #if necessary, add garamond font
fonts()

ggplot(df, aes(x= included, ymin =ymin, lower = low, middle = mid, upper = top, ymax = ymax, color = factor(included)))+
  geom_boxplot(stat = "identity")+
  theme_bw()+
  geom_hline(yintercept = 0, color= "black", lty= 2)+
  xlab("Inclusion at the NSS")+
  ylab("Discussion of the NSS the NSS")+
  ggtitle("Effect of Inclusion on Discussion of the NSS at the IAEA")+
  scale_color_manual("Status at NSS", values = c("violetred3","skyblue2"), labels=c("Excluded","Included"))+  theme(plot.title = element_text(family="EB Garamond", size=12),                                    legend.text = element_text(family="EB Garamond", size=10),legend.title = element_text(family="EB Garamond", size=12),axis.title = element_text(family="EB Garamond", size=12), axis.text = element_text(family="EB Garamond", size=10))

ggsave("output/rhetoric_boxplots.png") # creates figure 2

```


#Section 4: Matching Global Norms
Includes figures for appendix A2.2, A2.4, and generates data for Figure 5 and Table 5
```{r}
#_________________________________________
# Main model
mout1 = matchit(NSS2012.2016.x~ NWS+NPP+idealptavg+ GDP+ NAM + CTBTrat+ PSI + GTRI+ GICNT+ GP+ ICSANT+ CPPNM2005+ NWFZ +NATO+ NumSupplyGp + Region.y_1 + Region.y_2 + Region.y_3 + Region.y_4 + Region.y_5 + Region.y_6 + Region.y_7 + Region.y_8,  data = nss3, method= "genetic", pop.size = 5000, replace = TRUE)

bal.tab(mout1)# balance table, in Appendix A2.2

#running regressions
# With interaction
z.out_base = zelig(globalnorms~ NSS2012.2016.x, data = match.data(mout1), model = "ls", robust=T, cluster="ccode")

z.out = zelig(globalnorms~ NSS2012.2016.x*year + NWS+ NPP+ log(GDP) + NAM + CTBTrat+ PSI + GTRI+ GICNT+ GP+ ICSANT+ CPPNM2005+ NWFZ +NATO+ NumSupplyGp +  Region.y_1 + Region.y_2 + Region.y_3 + Region.y_4 + Region.y_5 + Region.y_6 + Region.y_7+USdist, data = match.data(mout1), model = "ls", robust=T, cluster="ccode")

# Without interaction
z.outB = zelig(globalnorms~ NSS2012.2016.x+ year + NWS+ NPP+ log(GDP) + NAM + CTBTrat+ PSI + GTRI+ GICNT+ GP+ ICSANT+ CPPNM2005+ NWFZ +NATO+ NumSupplyGp +  Region.y_1 + Region.y_2 + Region.y_3 + Region.y_4 + Region.y_5 + Region.y_6 + Region.y_7+ +USdist  , data = match.data(mout1), model = "ls", robust=T, cluster="ccode")
```

#Primary Stargazer Models (Table 5)
```{r}
stargazer( lm(z.out_base$zelig.out$z.out[[1]]), lm(z.outB$zelig.out$z.out[[1]]),lm(z.out$zelig.out$z.out[[1]]),
           dep.var.labels=c("Global Norms Score","Global Norms Score"),
           covariate.labels=c("Included at NSS","Year","Nuclear Weapons State","Nuclear Power","UN Ideal Avg","GDP","Non-aligned Movement","CTBT ratification","PSI","GTRI","GICNT","ICSANT","CPPNM 2005 Amendment","NWFZ","NATO","Number of Supply Groups","Africa","Asia","Europe-EU","Europe-NonEU","Latin America","North America","MENA","Oceana", "Included at NSS:Year"), align=TRUE,title= "Zelig Model Regression Results", nospace= TRUE , single.row = TRUE, star.cutoffs = c(0.05, 0.01, 0.001), notes = c(" * p<0.05; ** p<0.01; *** p<0.001"), notes.append = F, digits = 2)

```

##Figure 5 (right) Global Norms Boxplot
```{r}

########__________Figure 5 (right) raw data
# Simulating Values of global norms for included and excluded states, with year interaction.
#There is no neat way to transfer these, so I've manually input them below
x.out <- setx(z.out, NSS2012.2016.x=0)
x1.out <- setx(z.out, NSS2012.2016.x=1)
s.out <- sim(z.out, x = x.out, x1 = x1.out)
summary(s.out) #manually imported into Figure 3A, with code below

# Appendix plot A2.4- figure 1
pdf(file="Output/zeligFigure1.pdf")
plot(s.out)# Appendix plot
dev.off()



##########________ Appendix Plot A2.4- figure 2
#Simulated Results without interaction of year and inclusion
x.outB <- setx(z.outB, NSS2012.2016.x=0)
x1.outB <- setx(z.outB, NSS2012.2016.x=1)
s.outB <- sim(z.outB, x = x.outB, x1 = x1.outB)
summary(s.outB)

pdf(file="Output/zeligFigure2.pdf")
plot(s.outB)# Appendix plot
dev.off()


# manually exporting the zelig info for a box and whisker plot
mat = matrix(NA, 2, 2)


# expected values
mat[,1]= c(54.09114, 62.9109) # These values pulled directly from the s.out and s.outB summaries
#Standard Deviation
mat[,2] = c(1.61899, 0.6629443)

min(nss3$globalnorms[nss$NSS2012.2016.x==0], na.rm=TRUE)# min unincluded
max(nss3$globalnorms[nss$NSS2012.2016.x==0], na.rm=TRUE) # max unincluded

min(nss3$globalnorms[nss$NSS2012.2016.x==1], na.rm=TRUE)# min included
max(nss3$globalnorms[nss$NSS2012.2016.x==1], na.rm=TRUE) # max included

# These values pulled directly from the s.out summaries
df = as.data.frame(mat)
df$ymin = c(0,6)
df$ymax = c(100,100)
df$low = c(51.37401, 62.03917) # 2.5 %
df$mid = c(54.41826, 63.31403) # 50%
df$top = c(57.71755, 64.48857) #97.5 %
df$included = c(0,1)

ggplot(df, aes(x= factor(included), ymin =ymin, lower = low, middle = mid, upper = top, ymax = ymax, color = factor(included)))+
  geom_boxplot(stat = "identity")+
  theme_bw()+
  geom_hline(yintercept = 0, color= "black", lty= 2)+
  xlab("Inclusion at the NSS")+
  ylab("NTI Global Nuclear Security Norms Score")+
  ggtitle("Effect of Inclusion on Global Nuclear Security Norms")+
  scale_color_manual("Status at NSS", values = c("violetred3","skyblue2"), labels=c("Excluded","Included"))+
  scale_x_discrete(labels= c("Excluded", "Included"))+
  theme(plot.title = element_text(family="EB Garamond", size=12),                                    legend.text = element_text(family="EB Garamond", size=10),legend.title = element_text(family="EB Garamond", size=12),axis.title = element_text(family="EB Garamond", size=12), axis.text = element_text(family="EB Garamond", size=10))

ggsave("Output/norms_boxplots.png")

```

## Figure 5 (left) Global Norms through time graph 
Includes Figre 5B and Appendix figure A2.6
```{r}
# Overall, included states improve by much more (~30 pts as compared to ~8pts)

nss_avgs = nss3 %>% group_by(NSS2012.2016.x, year) %>% summarize(meannorm = mean(globalnorms, na.rm=TRUE))

# Resetting the year to make it easier to read in the graph
nss_avgs$year2 = nss_avgs$year+2011
nss3$year2 = nss3$year+2011


ggplot(nss_avgs, aes(x= year2, y = meannorm, color = factor(NSS2012.2016.x)))+
  geom_jitter(data = nss3, aes(x= year2, y = globalnorms, color = factor(NSS2012.2016.x)), alpha = 0.2)+
  geom_line()+
  theme_bw()+
  ylab("Mean Global Norm Scores")+
  xlab("Year")+
  scale_color_manual(name = "NSS", labels= c("Excluded","Included"), values = c("indianred3","lightskyblue"))+
  theme(plot.title = element_text(family="EB Garamond", size=12),                                    legend.text = element_text(family="EB Garamond", size=10),legend.title = element_text(family="EB Garamond", size=12),axis.title = element_text(family="EB Garamond", size=12), axis.text = element_text(family="EB Garamond", size=10))
  
ggsave("Output/globlnorms_time_all.png")

# ALLIES ONLY (appendix figure A2.6)

# However, there is selection in who is included

# looking only at NW Allies
allies_avgs = nss3[nss3$NWS==1,] %>% group_by(NSS2012.2016.x, year) %>% summarize(meannorm = mean(globalnorms, na.rm=TRUE))

# Resetting the year to make it easier to read in the graph
allies_avgs$year2 = allies_avgs$year+2011

ggplot(allies_avgs, aes(x= year2, y = meannorm, color = factor(NSS2012.2016.x)))+
  geom_jitter(data = nss3[nss3$NWS==1,], aes(x= year2, y = globalnorms, color = factor(NSS2012.2016.x)), alpha = 0.2)+
  geom_line()+
  theme_bw()+
  ylab("Mean Global Norm Scores")+
  xlab("Year")+
  scale_color_manual(name = "NSS", labels= c("Excluded","Included"), values = c("indianred3","lightskyblue"))+
  theme(plot.title = element_text(family="EB Garamond", size=12),                                    legend.text = element_text(family="EB Garamond", size=10),legend.title = element_text(family="EB Garamond", size=12),axis.title = element_text(family="EB Garamond", size=12), axis.text = element_text(family="EB Garamond", size=10))

ggsave("Output/globlnorms_time_allies.png")  
  
```


##Appendix A2.3: Robustness to a variety of matching algorithms 
```{r}
#____________________Additional Matching algorithms: for robustness ____________________________________

#trying all the different matching methods, both for balance statistics and to replicate the main result.

mout1_replace2 = matchit(NSS2012.2016.x~ NWS+NPP+IdealPoint +idealptavg+ GDP+ NAM + CTBTrat+ PSI + GTRI+ GICNT+ GP+ ICSANT+ CPPNM2005+ NWFZ +NATO+ NumSupplyGp + Region.y_1 + Region.y_2 + Region.y_3 + Region.y_4 + Region.y_5 + Region.y_6 + Region.y_7 + Region.y_8,  data = nss3, method= "genetic", pop.size = 5000, replace = TRUE, ratio = 2)

mout1_nearest = matchit(NSS2012.2016.x~ NWS+NPP+IdealPoint +idealptavg+ GDP+ NAM + CTBTrat+ PSI + GTRI+ GICNT+ GP+ ICSANT+ CPPNM2005+ NWFZ +NATO+ NumSupplyGp + Region.y_1 + Region.y_2 + Region.y_3 + Region.y_4 + Region.y_5 + Region.y_6 + Region.y_7 + Region.y_8,  data = nss3, method= "nearest",  replace = TRUE, ratio = 2)

mout1_optimal = matchit(NSS2012.2016.x~ NWS+NPP+IdealPoint +idealptavg+ GDP+ NAM + CTBTrat+ PSI + GTRI+ GICNT+ GP+ ICSANT+ CPPNM2005+ NWFZ +NATO+ NumSupplyGp + Region.y_1 + Region.y_2 + Region.y_3 + Region.y_4 + Region.y_5 + Region.y_6 + Region.y_7 + Region.y_8,  data = nss3, method= "optimal")


z.out_nearest = zelig(globalnorms~ year + NSS2012.2016.x + NWS+NPP+IdealPoint +idealptavg+ GDP+ NAM + CTBTrat+ PSI + GTRI+ GICNT+ GP+ ICSANT+ CPPNM2005+ NWFZ +NATO+ NumSupplyGp + Region.y_1 + Region.y_2 + Region.y_3 + Region.y_4 + Region.y_5 + Region.y_6 + Region.y_7, data = match.data(mout1_nearest), model = "ls")

z.out_replace2 = zelig(globalnorms~ year + NSS2012.2016.x + NWS+NPP+IdealPoint +idealptavg+ GDP+ NAM + CTBTrat+ PSI + GTRI+ GICNT+ GP+ ICSANT+ CPPNM2005+ NWFZ +NATO+ NumSupplyGp + Region.y_1 + Region.y_2 + Region.y_3 + Region.y_4 + Region.y_5 + Region.y_6 + Region.y_7, data = match.data(mout1_replace2), model = "ls")

z.out_optimal = zelig(globalnorms~ year + NSS2012.2016.x + NWS+NPP+IdealPoint +idealptavg+ GDP+ NAM + CTBTrat+ PSI + GTRI+ GICNT+ GP+ ICSANT+ CPPNM2005+ NWFZ +NATO+ NumSupplyGp + Region.y_1 + Region.y_2 + Region.y_3 + Region.y_4 + Region.y_5 + Region.y_6 + Region.y_7, data = match.data(mout1_optimal), model = "ls")

```

### Alternative Matching Algorithms - Appendix table
```{r}
#zout nearest #zout optimal #zout replace 2
stargazer( lm(z.out_nearest$zelig.out$z.out[[1]]), lm(z.out_optimal$zelig.out$z.out[[1]]),lm(z.out_replace2$zelig.out$z.out[[1]]),
           dep.var.labels=c("Nearest Neighbot","Optimal Matching","Genetic, 2 matches"),
           covariate.labels=c("Included at NSS","Year","Nuclear Weapons State","Nuclear Power","UN Ideal Avg","GDP","Non-aligned Movement","CTBT ratification","PSI","GTRI","GICNT","ICSANT","CPPNM 2005 Amendment","NWFZ","NATO","Number of Supply Groups","Africa","Asia","Europe-EU","Europe-NonEU","Latin America","North America","MENA","Oceana", "Included at NSS:Year"), align=TRUE,title= "Zelig Model Regression Results", nospace= TRUE , single.row = TRUE, star.cutoffs = c(0.05, 0.01, 0.001), notes = c(" * p<0.05; ** p<0.01; *** p<0.001"), notes.append = F, digits = 2)

```


## A2.8 Robustness Check: Replication in 2014-2020 controlling for 2012 Global Norms Score
```{r}

mout_later = matchit(NSS2012.2016.x~ NWS+NPP+idealptavg+ GDP+ NAM + CTBTrat+ PSI + GTRI+ GICNT+ GP+ ICSANT+ CPPNM2005+ NWFZ +NATO+ NumSupplyGp + Region.y_1 + Region.y_2 + Region.y_3 + Region.y_4 + Region.y_5 + Region.y_6 + Region.y_7 + Region.y_8,  data = nss3[nss3$year>1,], method= "genetic", pop.size = 5000, replace = TRUE)

z.out_2012control = zelig(globalnorms~ NSS2012.2016.x*year + NWS+ NPP+ log(GDP) + NAM + CTBTrat+ PSI + GTRI+ GICNT+ GP+ ICSANT+ CPPNM2005+ NWFZ +NATO+ NumSupplyGp +  Region.y_1 + Region.y_2 + Region.y_3 + Region.y_4 + Region.y_5 + Region.y_6 + Region.y_7  +gn2012+USdist, data = match.data(mout_later), model = "ls", robust=T, cluster="ccode")

z.outB_2012control = zelig(globalnorms~ NSS2012.2016.x+ year + NWS+ NPP+ log(GDP) + NAM + CTBTrat+ PSI + GTRI+ GICNT+ GP+ ICSANT+ CPPNM2005+ NWFZ +NATO+ NumSupplyGp +  Region.y_1 + Region.y_2 + Region.y_3 + Region.y_4 + Region.y_5 + Region.y_6 + Region.y_7  +gn2012+USdist, data = match.data(mout_later), model = "ls", robust=T, cluster="ccode")


# Stargazer output for A2.8 Regression table
stargazer(lm(z.outB_2012control$zelig.out$z.out[[1]]), lm(z.out_2012control$zelig.out$z.out[[1]]),  dep.var.labels=c("Global Norms Score"), align=TRUE, title= "Zelig Model Regression Results", nospace= TRUE , single.row = TRUE, star.cutoffs = c(0.05, 0.01, 0.001), notes = c(" * p<0.05; ** p<0.01; *** p<0.001"), notes.append = F, digits = 2)

```







